Load packages required for analysis
Load the fitted model and prior draws.
if (class(model) %in% "character") {
model <- ModelTBBCGEngland::read_libbi(file.path(model_dir, "libbi", "posterior"))
}
priors <- readRDS(file.path(model_dir, "data", "prior-params.rds"))
model
## Wrapper around LibBi
## ======================
## Model: Baseline
## Run time: 2070.558 seconds
## Number of samples: 1000
## State trajectories recorded: E H L P S T_E T_P YearlyAgeCases YearlyCases YearlyEAgeCases YearlyECases YearlyEPulCases YearlyEPulDeaths YearlyNonUKborn YearlyPAgeCases YearlyPCases YearlyPulCases YearlyPulDeaths
## Observation trajectories recorded: YearlyAgeInc YearlyHistPInc YearlyInc
## Parameters recorded: alpha_t_decay alpha_t_init c_eff c_hist chi_init delta epsilon_h_0_4 epsilon_h_15_89 epsilon_h_5_14 epsilon_l_0_4 epsilon_l_15_89 epsilon_l_5_14 HistMeasError kappa_0_4 kappa_15_89 kappa_5_14 M mu_e_0_14 mu_e_15_59 mu_e_60_89 mu_p_0_14 mu_p_15_59 mu_p_60_89 nu_e_0_14 nu_e_15_89 nu_p_0_14 nu_p_15_89 phi_0_14 phi_15_59 phi_60_89 rho_0_14 rho_15_59 rho_60_89 Upsilon_0_14 Upsilon_15_59 Upsilon_60_89 zeta_e_0_14 zeta_e_15_59 zeta_e_60_89 zeta_p_0_14 zeta_p_15_59 zeta_p_60_89
traces <- coda::mcmc(rbi::get_traces(model))
coda::rejectionRate(traces)
## M c_eff c_hist delta
## 0.992993 0.992993 0.992993 0.992993
## epsilon_h_0_4 epsilon_h_5_14 epsilon_h_15_89 kappa_0_4
## 0.992993 0.992993 0.992993 0.992993
## kappa_5_14 kappa_15_89 epsilon_l_0_4 epsilon_l_5_14
## 0.992993 0.992993 0.992993 0.992993
## epsilon_l_15_89 phi_0_14 phi_15_59 phi_60_89
## 0.992993 0.992993 0.992993 0.992993
## Upsilon_0_14 Upsilon_15_59 Upsilon_60_89 rho_0_14
## 0.992993 0.992993 0.992993 0.992993
## rho_15_59 rho_60_89 nu_p_0_14 nu_p_15_89
## 0.992993 0.992993 0.992993 0.992993
## nu_e_0_14 nu_e_15_89 zeta_p_0_14 zeta_p_15_59
## 0.992993 0.992993 0.992993 0.992993
## zeta_p_60_89 zeta_e_0_14 zeta_e_15_59 zeta_e_60_89
## 0.992993 0.992993 0.992993 0.992993
## mu_p_0_14 mu_p_15_59 mu_p_60_89 mu_e_0_14
## 0.992993 0.992993 0.992993 0.992993
## mu_e_15_59 mu_e_60_89 chi_init alpha_t_init
## 0.992993 0.992993 0.992993 0.000000
## alpha_t_decay HistMeasError
## 0.992993 0.992993
coda::autocorr.plot(traces)
- Evaluate Posterior Traces
plot_param(model, scales = "free", plot_type = "trace", burn_in = 0) + theme(legend.position = "none")
plot_param(model, prior_params = priors, scales = "free")
param_sum <- plot_param(model, prior_params = priors, plot_data = FALSE) %>%
group_by(Distribution, parameter, length) %>%
summarise(median = median(value),
lll = quantile(value, 0.025),
hhh = quantile(value, 0.975)) %>%
mutate(value = pretty_ci(median, lll, hhh, sep = ", ")) %>%
group_by(Distribution, parameter) %>%
mutate(Parameter = case_when(max(length) > 1 ~ paste(parameter, 1:n(), sep = "-"),
TRUE ~ as.character(parameter))
) %>%
ungroup %>%
select(Distribution, Parameter, value) %>%
spread(key = "Distribution", value = "value") %>%
select(Parameter, Prior, Posterior)
saveRDS(param_sum, file.path(model_dir, "data", "sum_prior_posterior.rds"))
knitr::kable(param_sum)
| Parameter | Prior | Posterior |
|---|---|---|
| alpha_t_decay | -0.13 (-0.23, -0.04) | -0.07 (-0.13, -0.04) |
| alpha_t_init | 0.84 (0.76, 0.90) | 0.79 (0.13, 0.87) |
| c_eff | 2.42 (0.13, 4.87) | 2.21 (0.16, 4.44) |
| c_hist | 12.57 (10.14, 14.82) | 12.99 (10.71, 13.16) |
| chi_init | 0.19 (0.08, 0.30) | 0.14 (0.08, 0.25) |
| delta | 0.78 (0.70, 0.86) | 0.78 (0.73, 0.80) |
| epsilon_h_0_4 | 1.53 (0.08, 4.67) | 0.49 (0.20, 5.01) |
| epsilon_h_15_89 | 23.87 (1.15, 74.94) | 20.29 (19.09, 57.16) |
| epsilon_h_5_14 | 3.66 (0.23, 11.82) | 3.50 (0.37, 14.25) |
| epsilon_l_0_4 | 581.14 (39.24, 1657.56) | 523.98 (191.49, 802.66) |
| epsilon_l_15_89 | 1106.77 (53.74, 3227.33) | 98.68 (78.38, 101.18) |
| epsilon_l_5_14 | 522.01 (29.10, 1523.41) | 311.46 (149.90, 1227.31) |
| HistMeasError | 1.00 (0.63, 1.42) | 1.16 (0.89, 1.36) |
| kappa_0_4 | 0.80 (0.04, 2.49) | 1.13 (0.25, 1.85) |
| kappa_15_89 | 1.15 (0.07, 3.44) | 0.90 (0.12, 1.98) |
| kappa_5_14 | 0.98 (0.04, 3.09) | 0.34 (0.06, 1.25) |
| M | 0.25 (0.01, 0.48) | 0.27 (0.04, 0.41) |
| mu_e_0_14 | 276.37 (206.60, 336.17) | 242.32 (177.20, 268.74) |
| mu_e_15_59 | 188.30 (76.36, 317.88) | 237.15 (129.68, 311.23) |
| mu_e_60_89 | 35.34 (1.87, 102.13) | 2.22 (0.44, 10.87) |
| mu_p_0_14 | 244.82 (156.68, 328.74) | 242.42 (202.11, 261.82) |
| mu_p_15_59 | 85.76 (4.20, 250.49) | 100.64 (51.90, 139.53) |
| mu_p_60_89 | 51.88 (2.09, 165.39) | 72.92 (11.51, 103.42) |
| nu_e_0_14 | 0.17 (0.00, 1.40) | 0.00 (0.00, 1.54) |
| nu_e_15_89 | 0.32 (0.01, 1.83) | 0.33 (0.06, 1.15) |
| nu_p_0_14 | 0.11 (0.00, 0.68) | 0.06 (0.05, 0.19) |
| nu_p_15_89 | 0.25 (0.02, 1.10) | 0.07 (0.02, 1.10) |
| phi_0_14 | 0.58 (0.28, 1.01) | 0.57 (0.49, 0.75) |
| phi_15_59 | 0.62 (0.27, 1.19) | 0.42 (0.41, 0.80) |
| phi_60_89 | 0.61 (0.28, 1.09) | 0.79 (0.48, 0.86) |
| rho_0_14 | 0.30 (0.27, 0.34) | 0.30 (0.27, 0.32) |
| rho_15_59 | 0.65 (0.64, 0.66) | 0.66 (0.65, 0.66) |
| rho_60_89 | 0.54 (0.52, 0.55) | 0.54 (0.53, 0.55) |
| Upsilon_0_14 | 0.63 (0.61, 0.65) | 0.64 (0.61, 0.65) |
| Upsilon_15_59 | 0.71 (0.70, 0.71) | 0.70 (0.70, 0.71) |
| Upsilon_60_89 | 0.75 (0.74, 0.76) | 0.75 (0.74, 0.76) |
| zeta_e_0_14 | 123.90 (55.57, 191.79) | 135.23 (95.99, 154.38) |
| zeta_e_15_59 | 60.60 (3.08, 184.34) | 31.35 (19.77, 77.57) |
| zeta_e_60_89 | 133.27 (61.54, 208.20) | 102.43 (27.49, 153.23) |
| zeta_p_0_14 | 94.05 (16.33, 176.11) | 63.57 (42.05, 116.35) |
| zeta_p_15_59 | 82.12 (2.64, 255.22) | 10.10 (4.90, 206.49) |
| zeta_p_60_89 | 124.89 (20.35, 252.94) | 56.52 (38.38, 308.59) |
plot_state(model, summarise = TRUE)
plot_state(model, summarise = TRUE, start_time = 59)
plot_state(model, states = c("YearlyHistPInc", "YearlyInc"), summarise = TRUE, start_time = 49) -> p2
p2
pred_obs <- plot_state(model, states = c("YearlyHistPInc", "YearlyInc", "YearlyAgeInc"), summarise = FALSE, start_time = 59, plot_data = FALSE)
obs <- ModelTBBCGEngland::setup_model_obs() %>%
bind_rows(.id = "state") %>%
mutate(time = time + 1931)
pred_obs <- pred_obs %>%
dplyr::filter(Average == "median") %>%
mutate(pred_value = pretty_ci(Count, lll, hhh, digits = 0)) %>%
left_join(obs, by = c("state" = "state", "time" = "time", "age" = "age")) %>%
select(Observation = state, Age = age, Year = time, `Observed Incidence` = value, `Predicted Incidence` = pred_value) %>%
mutate(Year = Year)
sum_obs_table <- pred_obs %>%
dplyr::filter(Observation %in% c("YearlyInc", "YearlyHistPInc")) %>%
drop_na(`Observed Incidence`) %>%
select(-Age) %>%
mutate(Observation = case_when(Observation == "YearlyInc" ~ "All TB cases",
Observation == "YearlyHistPInc" ~ "Pulmonary TB cases"))
saveRDS(sum_obs_table, file.path(model_dir, "data", "sum_obs_table.rds"))
age_cases_table <- pred_obs %>%
dplyr::filter(Observation %in% c("YearlyAgeInc")) %>%
drop_na(`Observed Incidence`) %>%
select(-Observation) %>%
group_by(Year) %>%
mutate(Age = c(paste0(seq(0, 65, 5), "-", seq(4, 69, 5)), "70-89")) %>%
ungroup
saveRDS(age_cases_table, file.path(model_dir, "data", "age_cases_table.rds"))
knitr::kable(sum_obs_table)
| Observation | Year | Observed Incidence | Predicted Incidence |
|---|---|---|---|
| Pulmonary TB cases | 1990 | 3618 | 3448 (2734 to 3752) |
| Pulmonary TB cases | 1991 | 3596 | 3310 (2855 to 3585) |
| Pulmonary TB cases | 1992 | 3816 | 3297 (2816 to 3623) |
| Pulmonary TB cases | 1993 | 3961 | 3292 (2645 to 3569) |
| Pulmonary TB cases | 1994 | 3694 | 3167 (2690 to 3547) |
| Pulmonary TB cases | 1995 | 3711 | 3180 (2733 to 3431) |
| Pulmonary TB cases | 1996 | 3690 | 3048 (2616 to 3400) |
| Pulmonary TB cases | 1997 | 3947 | 3106 (2626 to 3449) |
| Pulmonary TB cases | 1998 | 3926 | 3091 (2527 to 3391) |
| Pulmonary TB cases | 1999 | 3827 | 2975 (2496 to 3211) |
| All TB cases | 2000 | 1803 | 2152 (2043 to 2521) |
| All TB cases | 2001 | 1866 | 1949 (1948 to 2503) |
| All TB cases | 2002 | 1833 | 2012 (1849 to 2403) |
| All TB cases | 2003 | 1685 | 1919 (1919 to 2297) |
| All TB cases | 2004 | 1776 | 1795 (1747 to 2201) |
knitr::kable(age_cases_table)
| Age | Year | Observed Incidence | Predicted Incidence |
|---|---|---|---|
| 0-4 | 2000 | 75 | 2 (0 to 3) |
| 5-9 | 2000 | 59 | 1 (0 to 8) |
| 10-14 | 2000 | 75 | 0 (0 to 13) |
| 15-19 | 2000 | 112 | 7 (3 to 12) |
| 20-24 | 2000 | 133 | 15 (15 to 21) |
| 25-29 | 2000 | 116 | 45 (32 to 49) |
| 30-34 | 2000 | 147 | 81 (64 to 86) |
| 35-39 | 2000 | 99 | 129 (108 to 163) |
| 40-44 | 2000 | 83 | 175 (146 to 228) |
| 45-49 | 2000 | 105 | 227 (214 to 265) |
| 50-54 | 2000 | 116 | 323 (253 to 329) |
| 55-59 | 2000 | 112 | 318 (318 to 382) |
| 60-64 | 2000 | 101 | 326 (309 to 423) |
| 65-69 | 2000 | 99 | 300 (297 to 397) |
| 70-89 | 2000 | 371 | 229 (187 to 251) |
| 0-4 | 2001 | 113 | 3 (1 to 5) |
| 5-9 | 2001 | 65 | 1 (0 to 10) |
| 10-14 | 2001 | 51 | 1 (0 to 10) |
| 15-19 | 2001 | 130 | 3 (2 to 11) |
| 20-24 | 2001 | 145 | 14 (9 to 17) |
| 25-29 | 2001 | 127 | 32 (27 to 44) |
| 30-34 | 2001 | 122 | 76 (54 to 77) |
| 35-39 | 2001 | 128 | 98 (84 to 122) |
| 40-44 | 2001 | 107 | 180 (153 to 198) |
| 45-49 | 2001 | 101 | 215 (209 to 266) |
| 50-54 | 2001 | 96 | 278 (246 to 356) |
| 55-59 | 2001 | 91 | 335 (310 to 363) |
| 60-64 | 2001 | 94 | 336 (306 to 379) |
| 65-69 | 2001 | 105 | 307 (300 to 363) |
| 70-89 | 2001 | 391 | 214 (184 to 262) |
| 0-4 | 2002 | 102 | 4 (1 to 4) |
| 5-9 | 2002 | 62 | 1 (1 to 13) |
| 10-14 | 2002 | 64 | 1 (0 to 19) |
| 15-19 | 2002 | 129 | 7 (3 to 7) |
| 20-24 | 2002 | 148 | 12 (10 to 24) |
| 25-29 | 2002 | 120 | 36 (27 to 47) |
| 30-34 | 2002 | 112 | 59 (54 to 69) |
| 35-39 | 2002 | 127 | 112 (109 to 155) |
| 40-44 | 2002 | 98 | 147 (141 to 185) |
| 45-49 | 2002 | 89 | 220 (200 to 254) |
| 50-54 | 2002 | 104 | 282 (235 to 344) |
| 55-59 | 2002 | 116 | 328 (278 to 381) |
| 60-64 | 2002 | 88 | 309 (286 to 389) |
| 65-69 | 2002 | 90 | 328 (273 to 390) |
| 70-89 | 2002 | 384 | 214 (199 to 262) |
| 0-4 | 2003 | 74 | 2 (2 to 5) |
| 5-9 | 2003 | 46 | 1 (0 to 19) |
| 10-14 | 2003 | 59 | 1 (1 to 13) |
| 15-19 | 2003 | 102 | 6 (4 to 9) |
| 20-24 | 2003 | 116 | 13 (6 to 17) |
| 25-29 | 2003 | 137 | 26 (25 to 39) |
| 30-34 | 2003 | 118 | 54 (54 to 73) |
| 35-39 | 2003 | 113 | 98 (82 to 120) |
| 40-44 | 2003 | 109 | 152 (128 to 180) |
| 45-49 | 2003 | 83 | 201 (181 to 267) |
| 50-54 | 2003 | 102 | 252 (246 to 305) |
| 55-59 | 2003 | 92 | 282 (276 to 364) |
| 60-64 | 2003 | 98 | 288 (281 to 357) |
| 65-69 | 2003 | 94 | 290 (279 to 348) |
| 70-89 | 2003 | 342 | 242 (199 to 298) |
| 0-4 | 2004 | 112 | 4 (0 to 4) |
| 5-9 | 2004 | 71 | 0 (0 to 17) |
| 10-14 | 2004 | 81 | 1 (0 to 10) |
| 15-19 | 2004 | 111 | 5 (2 to 6) |
| 20-24 | 2004 | 120 | 14 (13 to 16) |
| 25-29 | 2004 | 136 | 28 (17 to 32) |
| 30-34 | 2004 | 129 | 53 (35 to 62) |
| 35-39 | 2004 | 103 | 101 (67 to 126) |
| 40-44 | 2004 | 128 | 132 (131 to 143) |
| 45-49 | 2004 | 94 | 188 (187 to 242) |
| 50-54 | 2004 | 103 | 240 (223 to 280) |
| 55-59 | 2004 | 98 | 282 (281 to 342) |
| 60-64 | 2004 | 86 | 297 (277 to 385) |
| 65-69 | 2004 | 92 | 307 (279 to 358) |
| 70-89 | 2004 | 312 | 212 (211 to 280) |
plot_state(model, states = c("YearlyAgeInc"), start_time = 59) + theme(legend.position = "none")